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Abstract 



The Hankel transform has found a wide variety of apphcations related to 



A method for computing the Hankel transform is proposed whereby 
the letter is reduced to a sum by representing the integrand as a smooth 
function times a Bessel function. The smooth function is replaced by its 
wavelet decomposition with a basis such its scalar product with the Bessel 
function is calculated analytically. The result is represented as a series, 
with the coefficients strongly depending on the local behavior of the func- 

^^ , tion being transformed. The application of the method is demonstrated 

\ jr^ • by an example illustrated with plots. 

o' 

f— ^ . problems in mathematical physics that possess an axial symmetry (see [1]). 

■^ ' The transforms of the zeroth (n = 0) and first (n = 1) kind, 
-4— > 

g: Fn{p) ^ J f{r)Jn{pr)rdr, 

F. : oo" (1) 

> . f{r)^ J Fn[r)Jn{pr)pdp. 

X 

;^ ■ are used most frequently. This is explained by the fact that, in Heaviside's calcu- 

5t 1 lus, the Laplace operator as applied to a function is equivalent to multiplication 

in the Laplace transform space: 

^ + ;s) «-) " -"'^"(ri, {^^\i- ^) fM ^ -pv,(„. 

In practice, however, this method encounters a number of difficulties. One 
of them is that integrals of type (1) cannot frequently be evaluated analyti- 
cally in finite form. Moreover, on computer, the function being transformed is 
specified as a set of numbers. Thus, an important task is to develop numerical 
methods for computing the Hankel transform with allowance made for the fact 



that Jn {pr) is a rapidly oscillating function for large values of p (or r) . Two 
approaches to the solution of this problem are available in the literature. The 
first is the fast Hankel transform, which was originally proposed in [2] and has 
been improved up to the present (e.g., see Appendix B in [3] and the references 
therein). This approach involves making a substitution and scaling, which re- 
duce the problem to logarithmic coordinates and the fast Fourier transform. 
This method applies to tabulated functions, but it involves conventional errors 
arising when a nonperiodic function is replaced by its periodic extension. Addi- 
tionally, the method is sensitive to the smoothness of functions in the space of 
logarithmic coordinates. An efficient algorithm (also based on the fast Fourier 
transform) was proposed for commuting the discrete Hankel (Fourier-Bessel) 
transform in [4]. It was noted in [4] that the algorithm is applicable in the 
case of an integral transform approximated by a sum. Another method (see [5]) 
involves the representation of the integrand as the product of a smooth function 
approximated by a set of polynomials and an oscillating function such that the 
integral of its product with the polynomial can be computed analytically. The 
accuracy of this approach is limited primarily by the accuracy of the polynomial 
approximation of the function to be transformed. 

The basic idea behind our method is similar to that in [5]. However, in- 
stead of a single approximation of the smooth factor, we apply its multiscale 
decomposition with respect to a function basis that takes into account the local 
properties of the function being transformed. Moreover, the scalar product of 
basis functions must be calculated analytically. 

These conditions are satisfied by some spline wavelet decompositions, which 
have been intensively studied in recent years (e.g., see [6-8]). Let g{r) be a 
function in L'^{R). Then it can be represented as a series 

9{r) = J2Y1 djki'jkir), (2) 

jez kez 

where ^pjk (r) — 2^/^tp2^r — k are basis functions with step k and scale j that sat- 
isfy the conditions of spatial localization and self-similarity under scale changes. 
The latter means that all ipjk are formed by contracting and translating a single 
function ip. 

Since the direct and inverse Hankel transforms (1) arc symmetric, we restrict 
our consideration to the direct transform and set g{r) = f{r)r for simplicity. 
Substituting (2) into the integral gives 



Pn{p) = X! X! '^^'^ / i'jk{r)JnPrdr. (3) 

Consider piecewise linear wavelets. Suppose that the basis wavelet is the piece- 
wise linear function 

•r-^ r — mh r — (m + l)h 

V{r) = 2^ 7 Om+l 37 Om, 



where bm are the values at the junction points and h is the minimum step. 

The hmits of summation are determined by the type of the wavelet chosen. In 
particular, as a basis, one can use the Battle-Lemarie -0^'^ and Stromberg tp^*''^ 
orthogonal wavelets [6], which are an infinite set of linear segments, or one can 
apply semiorthogonal or biorthogonal bases with a compact support consisting 
of M segments. Examples of the latter are B-spline semiorthogonal wavelets [8], 
biorthogonal wavelets obtained by applying the lifting scheme [9], and BlaC- 
wavclcts (Blending of Linear and Constant) [10], which combine a piecewise 
linear and a piecewise constant function. Note that the Battle- Lemarie and 
Stromberg bases rapidly decay at infinity, even though they have no compact 
support. Thus, they can be truncated to M segments in practical computations. 
Applying the integration formulas of [11], we find that the Hankel transform is 
equivalent to a sum: 

Foip) = I E E ^'^^djk E ((&fc+i - h) [{k + m+ l).h{2-\k + m+ l)hp)- 

jez kez m 

-{k + m) J i{2^^ {k + m)ph)\ + [(fc + m + 1)5^ - {k + m)bk+i] x 
{(A; + TO + l)Jo(2-J(fc + TO + l)hp) - (fc + m) Jo(2-J(fc + ■m)ph) + 
I [(fc + TO + 1)D{2-^ {k + m + l)hp) - {k + m)D{2-^ [k + m)ph)\ }) 

(4) 

Flip) ^ I E E '^''^djk E ( "<'"°+,^''"°^ [(fc + TO + l)D{2-^ (fc + TO + l)hp)- 

jez kez m ^ 

-(fc + m)D{2^^{k + ■m)ph)\ + [(fc + to + l)6fe - (fc + m)bk+i\ x 
{(fc + TO + 1) Jo(2-J(fc + m)hp) -(fc + m)Joi2-^{k + to + l)ph)}) 

(5) 
where D(x) = Ho{x)Ji(x) — Hi(x)Jo{x) and -ffo.i arc Struvc functions [11]. 

Thus, if the function being transformed (which is completely characterized 
by its coefficients [6-8]) has a finite-form definite integral of its product with 
a linear function, then scries (4) and (5) are an exact representation of the 
Hankel transform as a sum. Note that for sufficiently smooth functions strongly 
varying only in small ranges of their arguments, most of the rf-coefficients in (4) 
and (5) are small for high levels of resolution j. This makes it possible to derive 
an approximate numerical Hankel transform from multiscale analysis, in which 
case compactly supported bases have a significant advantage, because they avoid 
additional errors associated with the finite length of the function ranges used 
in numerical computations. Suppose that the function to be transformed is 
specified as a set of A^ points, which are interpreted as its projection on an 
interval onto the space of scaling functions of resolution J < log2 N. Then (2) 
is replaced by the series 

J 
di"^) ^ ^^kdjktpjk(r) +^sjkipJk{r), 
j=o k 

where tpok is a scaling function of a coarser resolution. In the computations, one 
can use only those of the coefficients djk that are greater than a prescribed small 
number e. The accuracy of the approximation of 5 by a truncated function g^ 



[7] is defined as A = (\\g - 



II) ' =1 '^d*f.'ipjk I • Taking into account 

the orthogonality of the basis wavelets, we obtain A < euq , where uq is the 
number of discarded coefficients d% < s. The resulting set contains wide bands 
of zeros in the domain where the function is smooth and an increased number 
of coefficients in the domain where it undergoes abrupt changes. Since the basis 
wavelet functions are narrowly localized, this implies a high degree of adaptivity 
in choosing an interpolation grid. Series (4) and (5) preserve their form when 
the limits in j and k are changed or terms are added that correspond to scalar 
products with scaling functions representable in a similar form. 

Thus, the integral transforms have been reduced to sums. For small argu- 
ments of the Bessel function, they can be computed by conventional codes and, 
for large arguments, by applying algorithms developed, for example, in [4]. 

Numerical computations in the case where g{r) is defined analytically are 
interesting on their own. For this purpose, we can use an inverse version of the 
lifting scheme [9] for piecewise linear basis functions. Let a function be defined 
on the interval [0, 1] (this can always be achieved by suitable scaling). This 
interval is bisected (the crudest approximation) and the lifting scheme is applied 
to a set of three function values (at the endpoints and the middle of the interval) 
to obtain the first iterates 



700 =9 
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yo{x) = ^ \k(pok{x), VFo(x) = 7ooV'oo(a;), g{x) « Vq{x) + Wo{x). 



k=0 



If the detailed coefficient 700 is not small, then the set is supplemented with 
the values of g{r) calculated at the middle points of each half-interval, and 
the next iterates are computed as follows. Calculate the detailed coefficient of 
coarser level J: 

7jfe=.g[2-('^+'=)(2fc + l)j -Vo [2-(''+'=)(2fc + l)j -Wj-i [2-(^+'=)(2fc + l) 

By using this relation, the detailed coefficients of the preceding levels and 
the approximation coefficients are updated to give 



7jk = Ijk 



-Wj 2 



)-{J+i) 



(2fc + i; 



J = 0, J - 1, 



7jfc-7jfc-VK,7(2-('^+i)(2fc+l)) 



Here Wj means all corresponding levels of the wavelet transform of difference 
between g(r) and the transform of level J — 1. 

The current representation of the function being decomposed is formed by 



1 ,7 2^-1 

^Afc(pofc(a;), T4^j-(a;) =^ ^ 7jfei/'ifc(a;), gix) 

j=0 k=0 



Vo{x) 



k=0 



i=o 



Here, the scaling function and the basis wavelet are given by 

ifiokix) = A(x - fc), 0(x) = A(2.T - 1) - -A(a;) - -K{x - 1), 

where K{x) = max{0, 1 — |x|}. 

As an example, we consider a function that has an exact first-kind Hankel 
transform: 

m^re-\ ^i(p)-(^exp(-^ 

Let a — 0.4. Then the basic nonzero part of the function is concentrated on 
the interval [0, 8]. By making the substitution r* ~ r/8, this interval is trans- 
formed into a unit interval, and the above procedures for wavelet decomposition 
with updating and the Hankel transform (5) are applied to g{r*) = f{r*)r*. 

Figure 1 shows the plot of g{r*) (solid curve) and its wavelet representa- 
tions at the zeroth (dot-and-dash), first (dotted), and second (dashed) levels. 
The corresponding curves in Fig. 2 plot the exact and approximate Hankel 
transforms in these cases. It can be seen that even this small number of levels 
provides an adequate pattern (except for the zeroth approximation, which gives 
a rather crude representation of the original function because of using the func- 
tion value at the middle point of the interval) . The next level provides a good 
approximation for most of the curve, except for large values of p* . This domain 
is corrected by the second approximation. Figure 3 shows the absolute error for 
this level, A{p*), which is equal to the difference between the approximate (dot- 
ted) and exact (solid) solutions against the background of the exact transform. 
It can be seen that the absolute error is fairly small in the entire domain shown 
in the figure and does not increase sharply with p* . The computed solution is as 
accurate as that obtained by linear approximation at the same nodes. However, 
our method provides rough estimation over the entire domain of the function, 
with subsequent refinement over intervals of interest. An improved accuracy in 
this approach can be achieved not only by increasing the number of decompo- 
sition levels taken into account but also by using a basis with unequally spaced 
nodes adapted to the behavior of the function to be transformed. 
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Figure 1: 
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Figure 2: 
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